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ABSTRACT 

To understand the basic mechanism of the formation of magnetic flux concentrations, we determine by 
direct numerical simulations the turbulence contributions to the mean magnetic pressure in a strongly stratified 
isothermal layer with large plasma beta, where a weak uniform horizontal mean magnetic field is applied. The 
negative contribution of turbulence to the effective mean magnetic pressure is determined for strongly stratified 
forced turbulence over a range of values of magnetic Reynolds and Prandtl numbers. Small-scale dynamo 
action is shown to reduce the negative effect of turbulence on the effective mean magnetic pressure. However, 
the turbulence coefficients describing the negative effective magnetic pressure phenomenon are found to be 
converged for magnetic Reynolds numbers between 60 and 600, which is the largest value considered here. In 
all these models the turbulent intensity is arranged to be nearly independent of height, so the kinetic energy 
density decreases with height due to the decrease in density. In a second series of numerical experiments, the 
turbulent intensity increases with height such that the turbulent kinetic energy density is nearly independent 
of height. Turbulent magnetic diffusivity and turbulent pumping velocity are determined with the test-field 
method for both cases. The vertical profile of the turbulent magnetic diffusivity is found to agree with what is 
expected based on simple mixing length expressions. Turbulent pumping is shown to be down the gradient of 
turbulent magnetic diffusivity, but it is twice as large as expected. Corresponding numerical mean-field models 
are used to show that a large-scale instability can occur in both cases, provided the degree of scale separation 
is large enough and hence the turbulent magnetic diffusivity small enough. 
Subject headings: MHD - Sun: magnetic fields - sunspots - turbulence 



1. INTRODUCTION 

In a stratified layer, magnetic fields do not normally 
stay in equilibrium but tend to become buoyantly unsta- 
ble (e.g . iNewcom bi 119611: IParkerl [196 6, 1979a; Gilman 
1970a b; Hughes & Proctor 1988; Cattaneo & Hughes 1988; 
Wissink et a l. 2000; Iso beet all 12005: Kersale et al. 2007), 
see also reviews by Hughes (2007) and Tobias & Weiss 
d2007l) . This mechanism related to magnetic buoyancy, is 
generally invoked in order to understand magnetic flux emer- 
gence at the solar surface (e.g. Hood et al. 2009). The mecha- 
nism does not explicitly rely upon the existence of turbulence, 
except that the origin of the Sun's m agnetic field is generally 
believed to be turbulent in nature; see Solanki et al] J2006I) for 
a recent review. 

Turbulent dynamos work under a variety of circumstances 
and are able to produce weakly nonuniform large-scale mag- 
netic fields (see Brandenburg & Subramanian 2005 for a re- 
view). At first glance this generation process is counter- 
intuitive, because it works against the well-kn own concept of 
turbulent mixing dTavlorill92lt iPrandtll [19251) . However, it 
is now well established that turbulence can also cause non- 
diffusive effects. In addition to the well-known a effect that 
is generally believed to be responsible for the Sun's large- 
scale field (IMoffattl [l97l IParkerl [19795 IKrause & R adler 
1980), there is also the A effect that is responsible for driv- 
ing the differential rotation o f the Sun jRudiger|[T980l 119891; 
Rudig er~& Hollerbachl 120041) . Yet another i mport effect is 
turbulent pumping or 7 effect (Radler 1969), which corre- 



sponds to the advection of mean magnetic field that is not 
associated with any material motion. The 7 effect appears, 
for example, in nonuniform turbulence and transports mean 
magnetic field down the gradient of turbulent intensity, which 
is usually downward in turbulent convection. However, this 
effe ct can also be modified by the mean magnetic field it - 
self dKitchatinov et al.H994rRogachevskii & Kleeorinl2006D . 
which can then correspond to a mean-field buoyancy effect. 

When invoking the concept of magnetic buoyancy, one 
must ask what the effect of turbulence is in this context. The 
turbulent pressure associated with the convective fluid mo- 
tions and magnetic fluctuations is certainly not negligible and 
reacts sensitively to changes in the background magnetic field. 
The main reason for this is that the kinetic energy density in 
isotropic turbulence contributes to the total turbulent dynamic 
pressure twice as mu c h as turbulent magnetic energy density 
dKleeorin et all [19891 [l990t iRogachevskii & Kleeorir] 120071 
hereafter referred to as RK07): 
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Here, Pturb is the total turbulent dynamic pressure caused by 
velocity and magnetic fluctuations, u and b, respectively, po 
is the vacuum permeability, p is the fluid density, and overbars 
indicate ensemble averaging. On the other hand, any rise in 
local turbulent magnetic energy density must be accompanied 
by an equal and opposite change of turbulent kinetic energy 



density in order to obey approximate energy conservation, i.e. 



\pu 2 



\b 2 /fi = E tot rj const. 



(2) 



Direct numerical sim ulations in open systems with bound- 
aries (Brandenburg et al. 2010, hereafter referred to as BKR) 
show that when the mean magnetic field B is much smaller 
than the equipartition field strength, B oq , the total energy is 
conserved, while when B < B cq , E tot decreases slightly 
with increasing mean field, and it varies certainly less than 
either pu 2 or b 2 ; see Fig. 1 of BKR. This clearly implies that, 
upon generation of magnetic fluctuations, the total turbulent 
dynamic pr essure shows a reversed (destabilizing) feedback 
jKleeorin et al.| [T990). i.e. 



Pturb = -^ 2 /^0 + 2S tot /3, 



(3) 



so both an increase of b 2 , as well as an increase of the imposed 
field, which decreases E to t, tend to lower the value of Pturb- 
For strongly anisotropic turbulence, Eq. ([3]) is also valid ex- 
cept for the change of the 1/6 factor into 1/2 (RK07). This 
phenomenology was supported by analytical studies using 
the spectral r relaxation approximation ( Kleeorin et al. 1990) 
and th e renormalization approach (Kleeorin & Rogachevskii 
1994) and led to the realization that the effective mean mag- 
netic pressure force (the sum of turbulent and non-turbulent 
contributions) is reduced and can be reversed for certain mean 
magnetic field strengths. Under certain conditions (e.g. strong 
density stratification), this can cause a magnetic buoyancy in- 
stability via perturba tions of a uniform me an magnetic field in 
stratified turbulence jKleeorin et al.ll 1 9931 1 1 996h . Later, when 
considering the effect of turbulent convection on the mean 
Lorentz force, RK07 suggested that magnetic flux concen- 
trations in the Sun such as active regions and even sunspots 
might be formed by this reversed feedback effect. 

Most of the numerical simulations on magnetic flux 
emergence (e.g.lStein & Nordj und 2001; Schussl er & Vogleri 



l2006riMartinezetal.ll2008t iRempel et all 120091) have been 
done using initial conditions with an already existing strongly 
inhomogeneous large -scale magnetic field. Recent simu- 
lations by BKR and IKitiashvili et all J2010I) ; iKapvlaetaT] 
(2012) study the formation of large-scale magnetic structures 
from an initially uniform large-scale magnetic field. In partic- 
ular, Large-Eddy Sim ulations of solar magneto-convection by 
IKitiashvili et al.l J201 0) give indications that the spontaneous 
formation of long-lived magnetic flux concentrations from an 
initial vertical uniform magnetic field might be possible, al- 
though the underlying mechanism in their simulations still re- 
mains to be clarified. A similar type of magnetic flux con- 
centration with vertical imposed field has b een seen incon- 
vection simulations at large aspect ratios bv lTao et~aI1 dl998l) . 
which show a segregation into magnetized and weakly mag- 
netized regions. One of the differences compared with BKR 
is the vertical orientation of the imposed magnetic field in 
turbulent convection. In forced a nd convection - driven tur- 
bulence simulations of BKR and Ka pyla et al.l J2012I) . re- 
spectively, the imposed magnetic field was a horizontal one. 
Other possibilities for causing flux concentrations include 
turbulent thermal collapse, whereby the magnetic field sup- 
presses the convective energy flux, leading to local cooling, 
and t hus to contraction and furth er enhancement of magnetic 
flux (iKitchatino v & Mazur 2000). By considering an isother- 
mal equation of state with isothermal stratification, we will 
exclude this possibility in our present work, allowing thus a 



more definitive identification of the effect of density-stratified 
turbulence on the mean Lorentz force. 

Meanwhile, direct numerical simulations of forced tur- 
bulence with an imposed horizontal magnetic field have 
demonstrated conclusively that in a simulation with an 
isothermal equation of state and an isothermal density 
stratification, spontaneous formation of mag netic structures 
dBrandenburg et al.l201 UlKemel et al.l20 12a) does indeed oc- 
cur. Those simulations used a scale separation ratio of 1:15 
and 1:30, i.e., the computational domain must be big enough 
to encompass at least 15 (or even 30) turbulent eddies in one 
coordinate direction. It does then become computationally 
expensive to achieve large Reynolds numbers, because their 
value is based on the size of the energy-carrying turbulent ed- 
dies rather than the size of the computational domain. In the 
present paper we restrict ourselves to a scale separation ratio 
of 1 :5 and are thereby able to demonstrate convergence of the 
turbulence coefficients describing the negative effective mag- 
netic pressure phenomenon for magnetic Reynolds numbers 
between 60 and 600. 

The basic phenomenon of magnetic flux concentration by 
the effect of turbulence on the mean Lorentz force has been 
studied by BKR based on numerical solutions of the mean- 
field momentum and induction equations. They demonstrated 
the existence of a linear instability for sufficiently strong strat- 
ification. This instability was followed by nonlinear satura- 
tion at near-equipartition strengths. Using direct numerical 
simulations (DNS) of forced turbulence, BKR also verified 
the validity of the phenomenology highlighted by Eq. (|3)- 
However, their DNS ignored the effects of stratification which 
would lead to additional effects such as turbulent pumping 
that might oppose the instability. 

Extending the DNS of BKR to the case with stratification is 
therefore one of the main goals of the present paper. This will 
allow us to make a meaningful comparison between DNS in 
a stratified fluid with mean-field modeling. We are now also 
able to present data for cases in which small-scale dynamo 
action is possible. This requires that the magnetic Reynolds 
number is large enough. As alluded to above, it is then ad- 
vantageous to choose a scale separation ratio that is not too 
extreme. While structure formation by the negative effective 
magnetic pressure phenomenon becomes impossible for small 
scale separation ratios, there is then also the advantage that 
the analysis becomes more straightforward in that horizon- 
tal and time averages can be employed. This would become 
problematic in the presence of structures that would break the 
assumptions of stationarity and homogeneity in the horizontal 
direction. 

There are two other possible caveats that may result from 
the simplification of using an isothermal equation of state with 
isothermal stratification. First, the effects of convection and a 
convectively unstable stratification on the mean Lorentz force 
are ignored. Fortunatel y, those turn out to b e weak, as shown 
in a separate paper by Kapyl a^et al.l (120121) . Second, owing 
to the spatial uniformity of the forcing function, u 2 is nearly 
uniform, so the effects of turbulent pumping down the gradi- 
ent of u 2 are ignored. We refer to these as models of type U. 
Furthermore, because of strong stratification of p, the equipar- 
tition field strength B cq = {papu 2 ) 1 / 2 also varies. This has 
the advantage that a single simulation with imposed field B 
spans a large range in the relevant control parameter Bo/B eq . 
However, in view of applications to turbulence in stellar con- 
vection this is unrealistic, because there pu 3 is nearly inde- 



pendent of height, so B eq increases with depth only like p 1 / 6 . 
For this reason we also study models in which B cq is nearly 
constant. We refer to these as models of type B. For this pur- 
pose we determine first the relevant turbulent pumping veloc- 
ity which is then used in a suitably adapted mean-field model. 
We begin by discussing first the determination of turbulent 
transport coefficients in Sect. [2] present our results in Sect. [3j 
focusing especially on models of type U, turn then in Sect. [4] 
to models of type B, compare the results at the level of mean- 
field models, and finish with a discussion of the main differ- 
ences between the magnetic buoyancy instability and the neg- 
ative effective magnetic pressure instability (NEMPI), before 
concluding in Sect. [5] 

2. DNS MODEL AND ANALYSIS 

We consider a cubic computational domain of size L 3 . The 
smallest wavenumber is then k\ = 2ir/L. We adopt an 
isothermal equation of state with constant sound speed c s , so 
the gas pressure is p = pc 2 s . The isothermal equation of state 
applies to both the background flow (the hydrostatic equilib- 
rium) and the fluctuating flow. In the presence of gravity, 
g = (0,0, —g), where g is the constant gravitational accel- 
eration, this leads to an exponentially stratified density, 



p= p exp{-z/H p ), 



(4) 



with a constant density scale height H p = c 2 /g and a nor- 
malization factor po- F or all our calculations we choose 
k\H p = 1. This implies that the number of scale heights 
is A In p = L/Hp = 2ir, corresponding to a density con- 
trast of exp 2tt m 535. This state is also chosen as our initial 
condition. Note that this is an equilibrium solution that is not 
affected by the possible addition of a uniform magnetic field 
B . 

We solve the equations of compressible magneto- 
hydrodynamics in the form 



DU 



= JxB-c 2 s Vp + V-(2vpS)+p(f + g), (5) 



dA 

— = UxB + V V 2 A, 
at 



dp 

at 



= -V-pU, 



(6) 



(7) 



where v and r\ are respectively kinematic viscosity and mag- 
netic diffusivity, B = Bo + V x A is the magnetic field con- 
sisting of an imposed uniform mean field, Bo — (0, Bo, 0), 
and a nonuniform part that is represented in terms of the mag- 
netic vector potential A, J = V x B/po is the current den- 
sity, and Sij = ^(diUj + djUi) — ^Sij'V ■ U is the traceless 
rate of strain tensor. The turbulence is driven with a forc- 
ing function / that consists of non-polarized random plane 
waves with an average wavenumber kf — 5 k\. The forc- 
ing strength is arranged such that the turbulent rms velocity, 
u r ms = {u 2 ) 1 / 2 , is around 0.1 c s . This value is small enough 
so that compressibility effects are confined to those associated 
with stratification alone. 

Our simulations are characterized by several non- 
dimensional parameters. We define the Reynolds number as 
Re = u vms /vkf, the magnetic Prandtl number as Vxm = ts/r) 
and the magnetic Reynolds number as Re^/ = RePrM- We 
anticipate that it is important to have Pr m < 1 ■ However, in 
order to reach somewhat larger values of Rejv/ we now choose 
as our primary model Pr^/ = 0.5 instead of 0.25, as was the 



case in BKR. In some additional cases, we span the entire 
range from Pim = 1/8 to Pr^/ = 8. For large enough values 
of Re m and Pim, there is small-scale dynamo action. We de- 
fine the equipartition field strength both as a function of z and 
for the middle of the domain, i.e. 



B eq (z) = (popu 2 ) ' , B cq0 = {papa) ' u 



1/2, 



(8) 



The latter will be used to specify the normalized strength of 
the imposed horizontal field, which is also independent of 
height. Another alternative is to normalize by the equipar- 
tition field strength at the top of the domain. In our models 
with nearly height-independent turbulent velocity, this would 
make the imposed field strength normalized by the equiparti- 
tion value at the top f=a 5 times bigger. 

In all cases we adopt stress-free perfect conductor boundary 
conditions at top and bottom of the domain. The simulations 
are performed with the PENCIL Code 1 , which uses sixth- 
order explicit finite differences in space and a third-order ac- 
curate time stepping method (Brandenburg & Dobler 2002). 

In this paper we present two groups of runs. In the first 
group we have the same forcing amplitude at all heights while 
in the second group we adjust the forcing such that the rms ve- 
locity depends on height such that the turbulent kinetic energy 
density is nearly independent of height. In contrast to earlier 
work where it was possible to analyze the results in terms of 
volume averages, we now have to restrict ourselves to hori- 
zontal averages which show a strong dependence on height. 
Thus, we determine the contribution to the mean momentum 
density that comes from the fluctuating field: 



£ , 

riy = ~puiuj + \kjb 2 - hbj, 



(9) 



where the po factor is dropped from now on and overbars indi- 
cate xy averages. The superscript f signifies the contributions 
from the fluctuating field. This, together with the contribution 
from the mean field, namely 



n 



pUiUj + S l3 (p+\B 2 ) - B t B - 2vpSij, (10) 



comprises the total mean momentum tensor, and the averaged 
momentum equation is given by: 



d_. 

di 



pUt = -v 



,(4 + < 



pgt- 



(ii) 



Here U and B are the mean velocity and magnetic fields, p is 
the mean fluid pressure. We are interested in the contribution 
to Eq. © that arises only from the presence of the mean field, 
so we subtract the corresponding tensor components that are 
obtained in the absence of the mean field. We thus define 



An' =rr M3 -n / ' 



(12) 



for which we make the following ansatz (RK07): 



An; 



■/ 



[\<h* 



q g gigj) B 2 + qsBiBj, (13) 



where g is the unit vector in the direction of gravity and the 
coefficients q p , q s and q g are expected to be functions of the 
modulus of the field, B = \B\. Equation ( TT3l can also be 
obtained from symmetry arguments, i.e., in the case of a hori- 
zontal imposed field, the linear combination of three indepen- 
dent true tensors, 6ij,gigj and BiBj, yields ansatz ( fT3l . 

' |http: //pencil-co de . googlecode ■ com| 



The meaning of the turbulence coefficients q p , q s and q g is 
as follows. The coefficient q p represents the isotropic turbu- 
lence contribution to the mean magnetic pressure, while q g 
is the anisotropic turbulence contribution to the mean mag- 
netic pressure, and the coefficient q s is the turbulence contri- 
bution to the mean magnetic tension. In the theory of RK07, 
the coefficients q p , q s and q g have been obtained using the 
spectral r approach and the renormalization approach. The r 
appro a ch has been justified in a number of numerical simula- 
tions ( Brandenburg et al. 2004J; iBrandenburg & Su bramanian 
2005b, 2007). However, if there is insufficient scale separa- 
tion, higher order terms such as Jjjj would need to be in- 
cluded. For helical flows, terms involving JjBj and BiJj, 
could also be present. Such terms are not included with the 
uniform fields used in the present study. 

The effective mean Lorentz force that takes into account the 
turbulence effects, reads: 



t= 300.0 




o 

-2 
-4 



-cM 



■pr?=-VMB*8 ii -B i B i 



« ■'•■Hj 



(14) 



Except for the contribution proportional to g^j and the fact 
that we use here horizontal averages, Eq. ( fl3l is equivalent 
to that used in BKR, where full volume averages were used. 
Asymptotic expressions for the B dependence of q p , q s , and 
q g are given in Appendix [B] Here we use DNS of density- 
stratified turbulence to determine these coefficients. In the 
present case, we have B sw (0, B, 0), so Eq. ( fT3l ) yields 



An4: 



-\q P B 2 , 



'{^q P -qs)B 2 , 



(15) 



All J zz = -(±q p + q g )B 2 
tt/ 



where we have computed AH U from DNS as 



ail{. 



P(u 2 



'o; 



) + I( 6 2_ 6 2 ) _ (6 2_ 6 2. )) 



(16) 



where no summation over the index i is assumed. The sub- 
scripts indicate values obtained from a reference run with 
B = 0. This expression takes into account small-scale 
dynamo action which can produce finite background mag- 
netic fluctuations 6o- (Thus, the reference run is not non- 
magnetic.) The critical magnetic Reynolds number for small- 
scale dynamo action is between 30 and 1 60, depending on the 
value of the magnetic Prandtl number (llskakov et al.l 120071: 
Brandenburg 201 1). Equation (TT3T > is then used to obtain ex- 
plicit expressions for 



q p = -2AU xx /B', 

q s = (Au f yy -Ani x )/W, 

q g = -(ATi f zz -ATll x )/B 2 , 



(17) 



allowing q p , q s , and q g to be evaluated at each height z. 



3. RESULTS 



3.1. 



Effective mean magnetic pressure 

We begin by considering the turbulence effects on the ef- 
fective mean magnetic pressure using a sequence of models 
of type U in which the rms velocity of the turbulence inten- 
sity is approximately independent of height, so B oq varies like 



FIG. 1. — Visualization of AB y /B cq on the periphery of the computa- 
tional domain for Model Ulh600 with B /B C q = 0.1 at Re M ~ 600 
using 512'^ meshpoints. 

p 1 / 2 and is about 23 times smaller at the top than at the bot- 
tom. In Fig. [T] we show a visualization of the departure of B y 
from the imposed field, AB y = B y — Bq, on the periphery of 
the computational domain for our model with the largest res- 
olution (Model Ulh600; for a complete list of all models dis- 
cussed in this paper see TableQ}. It turns out that most of the 
variability of the magnetic field occurs near the bottom of the 
computational domain. This is caused by the local variation 
of B eq ex p 1 / 2 . Therefore, Bq/ B oq is large in the upper parts, 
making it less easy for the turbulence to produce strong fluc- 
tuations due to the enhanced work done against the Lorentz 
force. By contrast, in the lower parts, Bo/B cq is small, al- 
lowing magnetic fluctuations to be produced. 

In the following we frequently use the symbol j3 to denote 
normalization by B cq , e.g., /3 = B/B cq (z). However, when 
we give the strength of the (z-independent) imposed or rms 
fields, we normalize with respect to B cq o, i.e., /3n. = Bo/B cq o 
and /3 Tms = B lms /B oq Q. The symbol /3 used here is not to 
be confused with the "plasma beta", which denotes the ratio 
of gas to magnetic pressures. To avoid confusion, we always 
spell out "plasma beta" in words. 

We have computed q p for all models of type U. We plot in 
Fig. |2] the dependence of q p on height for three different val- 
ues of Bo. In the following, the case with Bo/B cq o = 0.1 
will be used as our fiducial run. To improve the statistics, 
we present here time averaged results of q p , which itself is 
already averaged over x and y. Error bars have been calcu- 
lated by dividing the time series into three equally long pieces 
and computing the maximum departure from the total aver- 
age. In agreement with earlier work, q p is always positive 
and exceeds unity when the mean magnetic field is not suffi- 
ciently strong. This is the case primarily at the bottom of the 
domain (negative values of z) where the density is high and 
therefore the magnetic field, in units of the equipartition field 
strength, is weak. Since Bq = const and B cq increases with 
depth, Bo/B cq is smallest at the bottom, so q p also increases. 
The sharp uprise toward the lower boundary is just a result of 
the exponential increase of the density combined with the fact 
that the horizontal velocity reaches a local maximum on the 
boundary. 

The total effective magnetic pressure of the mean field 
(that takes into account the effects of turbulence on the mean 
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FIG. 2. — Dependence of q p (z) (solid lines) with error margins (dashed 
lines) as functions of z, for Models Ulh70, U2h70, and U5h70 with B = 
0.1B eq0 , Bo = 0.2B oq0 , and B = 0.5B cq0 , with Re M ss 70, g/c 2 a ki = 
1, and a density contrast of 530. Note that q v (z) reaches a maximum at the 
bottom of the domain where Bq / B cq (z) is minimal. 



Lorentz force) is given by ^[1 — q p (B)]B 2 . This has to be 

compared with the turbulent kinetic energy density, iptt 2 . 
Small contributions of terms oc q g to the effective mean mag- 
netic pressure are discussed in Sect. 3.3. In Fig.[3]we plot the 



TABLE 1 

Summary of all DNS models discussed in this paper. 



Model 



Re a/ Prj.j go q p o P+ 



Res. 



Ulh05 


0.5 


0.5 


0.1 


0.2 


0.32 


*0.04 


64 2 x 128 


Ulhl 


1.5 


0.5 


0.1 


0.2 


0.32 


*0.08 


64 2 x 128 


Ulh5 


5 


0.5 


0.1 


4 


0.34 


*0.15 


64 2 x 128 


UlhlO 


10 


0.5 


0.1 


13 


0.33 


*0.20 


128 3 


Ulh20 


23 


0.5 


0.1 


40 


0.33 


*0.29 


128 3 


Ulo35 


35 


1 


0.1 


90 


0.35 


*0.40 


128 3 


Ult35 


35 


2 


0.1 


70 


0.34 


*0.44 


128 3 


Ulf35 


35 


4 


0.1 


15 


0.28 


*0.44 


128 3 


Ule35 


35 


8 


0.1 


0.2 


0.32 


*0.38 


128 3 


Ulh40 


42 


0.5 


0.1 


170 


0.38 


*0.82 


128 3 


Ulq70 


70 


1/4 


0.1 


250 


0.37 


0.14 


128 3 


Ulh70 


70 


0.5 


0.1 


100 


0.33 


0.31 


128 3 


Ulh70h 


70 


0.5 


0.1 


60 


0.30 


0.29 


256 3 


Ulo70 


70 


1 


0.1 


50 


0.29 


0.42 


128 3 


Ult70 


70 


2 


0.1 


50 


0.26 


0.49 


128 3 


Ulf70 


70 


4 


0.1 


20 


0.22 


0.55 


128 3 


Ule70 


70 


8 


0.1 


- 


- 


0.49 


128 3 


U2h70 


70 


0.5 


0.2 


130 


0.35 


0.31 


128 3 


U5h70 


70 


0.5 


0.5 


200 


0.39 


0.31 


128 3 


Ulal40 


140 


1/8 


0.1 


200 


0.31 


0.43 


256 3 


Ulql40 


140 


1/4 


0.1 


40 


0.27 
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Here, /3o = Bq/B cci o and /3 rms = B rrnB /B ccl o denote field strengths in 
equipartition units, while /3* is a fit parameter that applies locally. Normally, 
/3rms = B TU1B /B cq refers to the field generated by small-scale dynamo ac- 
tion in the reference run with Bq = 0, except when there is an asterisk 
indicating that there is no small-scale dynamo and ft-ms gives the result from 
tangling of the applied field in the corresponding run with Bo ^ 0. Our 
reference model is indicated in bold face. 



FIG. 3. — Normalized effective mean magnetic pressure as a function of 
depth for the same models as in Fig. [2] Note that this function now reaches a 
negative minimum somewhere in the middle of the domain. 
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FIG. 4. — Same as Fig. [3] but as a parametric representation as func- 
tion of the local value of the ratio Bo/B C q(z). Note that the curves for 
Bq = 0.1B C qo, Bo = 0.2B cq o, and Bo = 0.5B cq o collapse onto a single 
dependency. The error range is the same as in the previous figure, but not 
shown for clarity. 



effective magnetic pressure normalized by B 2 



V cS = \(l-q p )B 2 /B, 



cq " 



(18) 



where B 2 itself is a function of height; see Eq. ®. It turns 
out that this function now reaches a negative minimum some- 
where in the middle of the domain. Work of Kernel et al. 
(2012b) has shown that the regions below the minimum value 
of "P c ff are those that can potentially display NEMPI. 

We expect that q p is a function of the ratio of B/B cq . This 
was observed numerically in BKR for constant B cq by vary- 
ing the value of B^ to obtain q p for a range of different sim- 
ulations. In the present case, however, B cq is a function of z, 
which is the main reason why q p depends on height. In Fig. [4] 
we plot the effective mean magnetic pressure as a function of 
magnetic field in units of the local equipartition value. Note 
that now all three curves for different values of B collapse 
onto a single curve, which demonstrates that the dependence 
of q p on both B and z can indeed be reduced to a single de- 
pendence on the ratio j3 = \B\/B oq (z). 

To quantify the form of the q p (f3) dependence, we used in 
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FIG. 5 . — Dependence of the fit parameters on Rsm using Models Ulh05- 
Ulh600. 



BKR a fit formula invo lving an arctan functi on. However, fol- 
lowing recent work of iKemel et al.l d2012bl) . a sufficient and 
certainly much simpler fit formula is 



QpW) 



P 2 IPl' 



(19) 



Here the fit parameters q p $ and f3 p are determined by mea- 
suring the minimum effective magnetic pressure, "P m i n = 
min('P ff), as well as the position of the minimum, i? m ; n , 
where V c s (-B m in) = 'Pmin- For our setups, 7> min is typically 



around —0.05, while /3 n 



Bmin/Beq is between 0.1 and 



0.2. This is remarkably close to Fig. 3 of RK07, who used 
the spectral r relaxation approximation. The fact that nearly 
the same functional form for the effective magnetic pressure 
of the mean field is obtained, supports the idea that this effect 
is robust. 

For many of the models in Table Q] we have determined the 
fit parameters q p o and f3 p . It turns out that for small values of 
Rej\/, q p0 increases quadratically with Re^ and f3 p decrease 

like Re^ 1 . Thus, for Re^/ < 30, we have ft 2 = /3 2 q p o ~ 
const. The significance of this is that /3* turns out to be nearly 
independent of Re^/; see Fig. [5] It allows rewriting the fit 



m 



pi 



Pi+p 2 



(20) 



where for small values of Re^/, only /3 p depends on Re^/. In 
particular, we have then j3 p ~ 1.05/Rejv/. 

The obtained scaling for q p o ~ O.lRe^j is consistent with 
an estimate based on the quasi -linear c alculations. This anal- 
ysis is similar to that of Rudig eret al.l J2012I) , except that we 
performed an explicit integration in wfc-space for a power-law 
kinetic energy spectrum of the background turbulence and for 
a Lorentz profile for the frequency dependence of the veloc- 
ity correlation function. For /? <C 1 this analysis yields the 
expression (see AppendixlAb 



Qpo 






60 Pr 



(21) 



M 



where Re^j/PrM = Rejv/Re. For Pr^/ = 1/2, we have q p — 
0.25 Re\.j, which is in qualitative agreement with our scaling 
for q P Q. The discrepancy in the coefficient is related to the fact 
that the quasi-linear approach is only valid for small magnetic 
and fluid Reynolds numbers. Therefore, the limit Pim — > in 
the framework of the quasi-linear approach only implies the 
case of large magnetic diffusion r\, while the case of small v 
need to be considered in the framework of approaches that are 
valid for large fluid Reynolds numbers (like the r relaxation 
approach). 

Looking at Table [T] it may seem surprising that q p o can 
reach values as large as 250 (see, e.g., Model Ulq70 with 
Rem ~ 250 and Ptm = 4). However, the more relevant 
quantity is q p (/3)/3 2 , which is of the order of 01 for large field 
strengths. Summarizing the results from Fig. [5] we find 



P P « 1.05 Re M \ = 0.33 (forRe M < 30), (22) 



and 



P P ra 0.035, 



,3* = 0.23 (for Re M > 60), 



(23) 



which implies that q p ((3)/3 2 is below 0.1 and 0.06 for the 
regimes applicable to Eqs. d22b and d23l . respectively. Also, 
it is tempting to associate the sudden drop of from 0.33 to 
0.23 with the onset of small-scale dynamo action. The fact 
that small-scale dynamo action reduces the negative effect of 
turbulence on the effective mean magnetic pressure was al- 
ready predicted by RK07 and is also quite evident by looking 
at Table [T] where q p o is found to reach more moderate values 
after having reached a peak at ReM ~ 30. 

We reiterate that, as long as the value of the plasma beta 
(i.e. the ratio of gas pressure to magnetic pressure), is much 
larger than unity, our results are independent of the plasma 
beta. What matters is the ratio of magnetic energy density to 
kinetic energy density, not the thermal energy density. This 
is also clear from the equations given in Appendix [B] In the 
present simulations, the plasma beta varies from between 5 
and 100 at the top to around 10 5 at the bottom, so the total 
pressure (gas plus magnetic plus turbulent pressure) is always 
positive. 

3.2. Dependence on magnetic Prandtl number 

In most of the runs discussed above we used Pr^/ = 1/2. 
As expected from earlier work of RK07, the negative mag- 
netic pressure effect should be most pronounced at small 
PrM- This is indeed confirmed by comparing with larger and 
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FIG. 6. — Normalized effective mean magnetic pressure for different values 
of Pr M , for Models Ulq70, Ulh70, Ult70, and Ule70, where Re M « 70 
and 73o = 0.1/3 GqU . 



FIG. 8. — Similar to Fig. [4] but for q s B 2 /B^ and now as a function of 
By /.B C q. Within the error range (dashed lines), q s = for all field strengths. 
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FIG. 7. — Resolution dependence of q p (B y /B C q) for Prjy = 1/2, 
Re M « 70 for Models Ulh70 and Ulh70h using 128 3 and 256 3 mesh 
points. Error bars are marked by the dashed lines. 

smaller values of Pi"m; see Fig. [6] where we show the results 
for Models Ulq70, Ulh70, Ult70, and Ule70. Here and in 
Table Q]the different values of Pi - a/ are denoted by the letters 
q, h, o, t, f, and e, which stand for a quarter, half, one, two, 
four, and eight, while the letter a is used for one eighth. 

3.3. Resolution dependence 

A density contrast of over 500 may seem rather large. How- 
ever, this impression may derive from experience with poly- 
tropic models (see, e.g.. lCattaneo et aiH l991). where most of 
the density variation occurs near the surface. In our isother- 
mal model, the scale height is constant, so the logarithmic 
density change is independent of height. Figure [7] shows that 
the error bars for the 256 3 run are smaller than those for the 
128 3 run, and that the minimum of V c s is somewhat less shal- 
low, but within error bars the two curves are still comparable. 
Here, the relevant input data are averaged over a time interval 
At such that u rms kfAt is at least 1500, and that error bars, 
which are based on 1/3 of that, cover thus at least 500 turnover 
times. 

3.4. Coefficients q s and q g 
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FIG. 9. — Similar to Fig. [8] but for q g . Note that q g is positive. 

Depending on the size and magnitude of the coefficients 
q s and q g , their effect on the instability could be signifi- 
cant. Most importantly, a positive value of q s was found to 
be chiefly respo nsible for producing three-dimensional mean- 
field structures (Kernel et al. 2012b), i.e., structures that break 
up in the direction of the imposed field. The coefficient q g , 
on the other hand, affects the negative effective magnetic 
pressure and could potentially enhance its effect significantly 
(RK07). 

Using Eqs. ( fT5l ). we now determine q s and q g . The results 
are shown in Figs.[8]and|9]for the three imposed field strengths 
considered above, where u 2 is nearly independent of z, so 
B eq (z) varies by a factor of exp7r w 23, allowing us to scan 
the dependence on B/B cq in a single run. It turns out that 
both q s and q g are around zero. 

Recent DNS of stratified convection with an imposed hori- 
zontal magnetic field did actually yield non-vanishing (pos- 
itive) values of q g for stratified convection (Kapvlaetal. 
12012b . In the present study with vertical density stratifica- 
tion, q g is much smaller, but generally positive. This ap- 
pears to be in conflict with the theoretical expectation for q g 
given in Appendix iBl where q g — 0(£ 2 f /H 2 ) if we assume 
If = 2-7r/fc/ w 1.3, which gives q g f=a 1.2. However, with- 
out the 2n factor, we would have if — kj 1 — 0.2, and thus 
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FIG. 10. — Visualization of AB y /B ccl on the periphery of the computational domain for runs with -Bo/^eqO = 0.1, 0.2, and 0.5, with Re = 70. 
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FIG. 1 1 . — Normalized mean magnetic field in the direction of the imposed 
field versus height for Bq = 0.1-B cq o, 0.2B cq o, and 0.5-B cq o with Re = 
70. 

q g s=s —0.03. Fig.[9]suggests a positive value of similar mag- 
nitude. This issue will hopefully be clarified soon in future 
work. 

Next, we discuss the results for q s . In BKR there was some 
evidence that q s can become positive in a narrow range of field 
strengths, but the error bars were rather large. The present re- 
sults are more accurate and suggest that q s B 2 /B 2 is essen- 
tially zero. This is also in agreement with recent convection 
simulations (Kapvla et al. 2012). 

In summary, the present simulations provide no evidence 
that the coefficients q s and q g could contribute to the large- 
scale instability that causes the magnetic flux concentrations. 
This is not borne out by the analytic results given in Ap- 
pendix [B] The results from recent convection simulations 
fall in between the analytic and numerical results mentioned 
above, because in those q g was found to be positive, while q s 
was still found to be small and negative (Kapyla e t al.ll2012h . 

4. COMPARISON OF MODELS OF TYPES U AND B 
4.1. Results from DNS 

As discussed in the beginning of Sect. lXTl most of the vari- 
ability of the magnetic field occurs near the bottom of the 
computational domain. This is also evident from Fig. [T0J 
where we show visualizations of the y component of the 
departure from the imposed field, AB y , on the periphery 
of the domain for Models Ulh70, U2h70, and U5h70 with 
Bo/B cq o — 0.1, 0.2, and 0.5, respectively. 

The vertical dependence of the horizontally averaged mean 
magnetic field (now normalized to Bq — const) is shown 



FIG. 12. — Turbulent kinetic energy density versus height for n = oo 
for Bo = (dotted line) and Bq = 0.2B oq o compared with the case for 
n = 1.4 and B = 0.2B cq0 . 



in Fig. QT| We see that, especially for intermediate field 
strengths, there is an increase of the magnetic field near the 
top of the domain. One possibility is that this is caused by 
the effect of nonlinear turbulent pumping, which might cause 
the mean field to be pumped up due to the gradients of the 
mean turbulent kinetic ener gy density in the presence of a fi- 
nite mean magnetic field (cf. lRogachevskii & Kleeorin 2006). 
This type of pumping is different from the regular pump- 
ing down the gradient of turbulent intensity (Radler 1969). 
To eliminate this effect, we have produced additional runs 
where the kinetic energy density is approximately constant 
with height. This is achieved by modulating the forcing func- 
tion by a z-dependent factor e z / H i . We define n = Hj/H p 
and find that for n = 1.4 the kinetic energy density is approx- 
imately independent of height; see Fig. [12] 

As a consequence of reducing the turbulent driving in the 
lower parts by having B cq (z) « const, we allow the magnetic 
field to have almost the same energy density as the turbulence, 
i.e. Bo/B cq (z) is approximately independent of z. This also 
means that the fluctuations are now no longer so pronounced 
at the bottom of the domain (Fig.[T3l>. where Re drops to val- 
ues around 5 and the flow is no longer turbulent. However, at 
the top the Reynolds number is around 120, so here the flow 
is still turbulent. In Fig.[l4]we show the vertical dependence 
of the horizontally averaged mean magnetic field in units of 
the imposed field strength. Note that now the field shows an 
increase toward the bottom of the domain. This effect might 
be related to regular turbulent pumping (Radler 1969), which 
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FIG. 13. — Visualization of AB y /B cq on the periphery of the computa- 
tional domain for Model Blh35 with nearly uniform turbulent kinetic energy 
density using Hi = nH p with n = 1.4. 
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FIG. 14. — Normalized mean magnetic field in the direction of the im- 
posed field versus height in the case of nearly constant turbulent kinetic en- 
ergy density, i.e. B Bq (z) ?s const, for Models B07h35 and B2h35 with 
Bo/B cq o = 0.07 and 0.2, respectively. 

now has a downward component because u 2 decreases with 
depth. 

4.2. Determination of j] t and 7 from the simulations 

We use the test-field method of ISchrinner et al.l (120051 
20071) in the Cartesian implementation, as described by 



Brandenb urg et al. l d2008al) . to compute 774 and 7 from the sim- 
ulations in the presence of the applied field. We refer to this as 
the quasi-kinematic test-field method, which is applicable if 
the magnetic fluctuations are just a consequence of the mean 
field; see Rheinhardt & Brandenburg (2010) for details and 
extensions to a fully nonlinear test-field method. For further 
comments regarding the test-field method see Appendix ICl 
We analyze the two setups discussed above, namely those of 
type U (where u rms and hence r\ t are nearly constant in height) 
and those of type B (where B cq is nearly constant in height). 
The set of test fields includes constant and linearly grow- 
ing ones. For both models we use Bq — O.Olp 1 / 2 ^, cor- 
responding to Bq m 0.1B cq o for models of type U and 
Bq sa 0.07-Bcqo for model of type B. The results are shown 
in Figs. G3] and [16] for Models B07h35 and Ulh70, respec- 



FIG. 15. — Results for r\t and 7 with the test-field method (solid lines; 
error margins are shown as dashed lines) for Model B07h35. In the upper 
panel, the dotted line gives 1.2u r ma/3fc^ and the dash-dotted line repre- 
sents 1.2 « ms i) exp(z/H u )/3kf. In the lower panel, the dash-dotted lines 
represents —1.2 « r ms0 exp(z/H u )/3kfH u . 
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FIG. 16.— Same as Fig.fTBI but for Model Ulh70. 

tively. In Table |2] we summarize the relevant parameters in- 
ferred for these models. It turns out that the DNS results are 
well described by r\ t = l-2?7to, with rjto = u lms /3kf and 
u Tms (z) = u rms0 exp(z/H u ), but 7 = -drjt/dz, i.e. without 
the 1/2 factor expected from the kinematic mean-field theory 
(lRadleil969h . 

We emphasize that 7 sa for models of type U, suggesting 
that additional effects owing to the mean magnetic field such 
as mean-field magnetic b uoyancy ( Kitchatinov et al. 1994; 
Roga chevskii & Kleeorinll2006l) are weak (AppendixlDt. 

4.3. Comparison with mean-field models 



10 



4.3.1. Basic equations 

We follow here the same procedure as BKR and consider 
the equations for the mean velocity U, the mean density p, 
and the mean vector potential A in the form 



dU 

~~dt 



-U-VU-c 2 s Vlnp + g + F 



dA 
~~dt 



M 



T 



KVtot) 



U x B 



ot 



- 7] J - V$, 

= -V-pU, 



(24) 



(25) 



(26) 



where $ is the mean electrostatic potential, B = Bq + V x A 
is the mean magnetic field including the imposed field, and 



■tPEM _ T w "D 



£F* =JxB + ±V(q p B z 



(27) 



is the effective mean Lorentz force, where we use for q p (B) 
the fit formula give n by Eq. ( fT9l ). However, in view of the 
results of Sect. 13.41 the q s and q g terms will now be omitted, 
and 



FiCtot = (vt + v) (V'C/ + VV • U + 2SV In 75) (28) 
is the total (turbulent and microscopic) viscous force, 



£ = 7 x B - TjtJ, 



(29) 



is the mean electromotive force, where 7 is the turbulent 
pumping velocity and rjt is the turbulent magnetic diffusivity. 
In ou r mean-field models we a ssume vt Ivt = 1 (JY ousef et al. 
120031) . The kinematic theory of lRoberts & Sowardld 19751) and 
others predicts that r]t(z) = u rms (z)/3kf and 7 = — 5 V774. 
It is fairly easy to assess the accuracy of these expressions by 
computing turbulent transport coefficients from the simula- 
tions using the test-field method ( Schrinner et al. 2005, 2007). 
A comment regarding $ is here in order. It is advantageous 
to isolate a diffusion operator of the form rjt V 2 A by using the 
so-called resistive gauge in which <J> = — 774 V- A. This means 
that the diffusion op erator now becomes 77 t V 2 A+(VA)Vrfe 
(Dobleret al. 2002). This formulation is advantageous in sit- 
uations in which r\ t is non-uniform. 

4.3.2. Results from the mean-field models 



Next we consider solutions of Eqs. (I24li-(|29jl for models of 
types U and B using the parameters specified in Table [2] To 
distinguish these mean-field models from the DNS results, we 
denote them by script letters U and B. We have either constant 
■qt (Models U) or constant B cq (Models B). In both cases we 
use )) = 2v = 4x 10~ 4 c s /fci, B = 0.005p 1/2 c s , with q^ = 
40 and j3 p = 0.07, while B cq o and po are given in Tableland 
correspond to values used in the DNS. This gives the profile 
of u lms (z) = B cq /yfp, which allows us to compute rj t (z) = 
u T m.s(z)/3kf and v t (z) = rjt(z) for an assumed value of kf. 
In the DNS presented here we used kf/k i = 5 and did not 
find evidence for NEMPI, but the DNS of iBrandenburg et all 
(201 1) and Kernel et al. (2012a) for kf/ki = 15 and kf/ki = 
30, respectively did show NEMPI, so we mainly consider the 
case kf/ki — 15, but we also consider kf/ki = 5 and 10. 

As in BKR, Eqs. (124-b — (f26b exhibit a linear instability with 
subsequent saturation. However, this result is still remarkable 
because there are a number of differences compared with the 
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FIG. 17. — Evolution of the mean velocity for Models hi and B obtained 
by solving the mean field equations. "Model £30 (7 = 0)" refers to a model 
where the pumping velocity is ignored. In all cases, kf/ki = 15 is assumed. 
In the inset we compare the evolution for Model hi with those for smaller 
values of fey/fei (Model U 10 and US). 



models studied in BKR. Firstly, we consider here an isother- 
mal atmosphere which is stably stratified, unlike the isen- 
tropic one used in BKR, which was only marginally stable. 
This underlines the robustness of this model and shows that 
this large-scale instability can be verified over a broad range 
of conditions. Secondly, this instability also works in situa- 
tions where rjt and/or B cq are non-uniform and where there 
is a pumping effect that sometimes might have a tendency to 
suppress the instability. 

In Fig. [17] we compare the evolution of the rms velocity of 
the mean flow. Note that, in contrast to the corresponding 
plots in BKR, we have here normalized U rms with respect to 
vao = Bo/yfpo and time is normalized with respect to the 
Alfven wave traveling time, (vAoki) -1 . This was done be- 
cause in these units the curves for Models U and B show sim- 
ilar growth rates. This is especially true when the pumping 
term is ignored in Model B0, where we have set artificially 
7 = 0. With pumping included (as was determined from 
the kinematic test-field method), the growth rate is slightly 
smaller (compare dashed and dotted lines). The pumping 
effect does not significantly affect the nonlinear saturation 
phase, i.e. the late-time saturation behavior for the two ver- 
sions of Model B is similar. Instead, the saturation phase is 
different for Model U compared with Model B and the satu- 
ration value is larger for Model U. The inset of Fig. [T7lcom- 
pares the results for Model U (with kf/ki — 15) with Models 
U10 and U5 for kf/ki = 10 and 5, respectively. Note that 
NEMPI is quite weak for kf/k\ = 5, which is consistent with 
the DNS presented here. 

Visualizations of the mean magnetic field as well as the 
mean velocity are shown in Fig. [18] for three different times 
near saturation for Model UQ2. Here, we use a weaker im- 
posed field, Bq = 0.002 / o 1 / 2 c s , corresponding to Bo/B eq p = 
0.02, so that NEMPI starts closer to the surface. For 
B /B c(i0 = 0.05, NEMPI starts in the middle of the do- 
main, leaving less space before the structures have reached 
the bottom of the domain. Note the converging flows toward 
the magnetic structures, with the largest velocities occurring 
in the upper layers where the density is smallest. 

4.3.3. Comments on the shape of mean-field structures 
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FIG . 18. — Mean magnetic field in the y direction (color coded) together 
with velocity vectors in the xz plane for Model W02. Note the spontaneous 
formation of flux structures. 

The descending structures found in the present mean-field 
calculations are qualitatively similar to those of BKR who 
considered a polytropic layer. In both cases the structures sink 
and become wider. This is quite different from the behavior 
of individual turbulent eddies and flux tubes that one tends to 
monitor in DNS. Clearly, individual magnetic structures ex- 
perience magnetic buoyancy, where the vertical motion is the 
result of a balance of magnetic buoyancy and downward ad- 
vection by the ambient flow (see Fig. 10 of Brandenburg et al. 
1996). The mean-field model cannot describe individual 
(small-scale) structures, although the net effect of the vertical 



TABLE 2 
Summary of parameters entering models U and B. 



Model 



'/(o 



B, 



cqO 



A) 



Re 



M 



comment 



U 

W10 

U5 

UQ2 

8 

BO 



15 
10 
5 

15 
15 
15 



0.0072 
0.0048 
0.0024 
0.0072 
0.0024 
0.0024 



0.11 

0.11 

0.11 

0.11 

0.036 

0.036 



0.05 
0.05 
0.05 
0.02 
0.05 
0.05 



36 

36 

36 

36 
2.5-60 
2.5-60 



B c . 



3cq0 I 



z/2H p 



B cq = B eq0 e-*/ 2H >> 



B c , 



VI 
7 = 



B oq I 
BcqO I 



z/2ff p 
-z/2ff„ 



Vtoe 



■ /2H P 



assumed 



Tildes indicate nondimensional quantities: 

1/2 
ritokl/Cs, B eq = Bcq0/P < 



kf/kx, fjto 



, while A) = B /B, 



cqO- 



motion of individual structures results in a turbulent pumping 
velocity 7 of mean-field structures. 

While turbulent downward pumping has been seen in nu- 
merous DNS, the negative effective magnetic pressure insta- 
bility is a new effect that has been seen so far only in the 
DNS o f forced turbulence in BKKMR; see also Kernel et al. 
d2012al) . Amazingly, those DNS show very similar structures 
resembling that of a descending "potato sack" (see Fig. 1 of 
BKKMR). They sink because the negative effective magnetic 
pressure is compensated by increased gas pressure, which in 
turn leads to larger density, so they become heavier than the 
surroundings. However, these turbulent magnetic structures 
are only poorly associated with material motion; see the flow 
vectors in Fig. [18] Therefore, a change of their volume is 
not governed by mass conservation. In particular, these struc- 
tures do not become narrower during their descent as individ- 
ual blobs do in a strongly stratified layer. 

We should point out that these mean-field structures do not 
always initiate at the top of the layer. The initiation height 
depends on the value of z where V e s reaches a minimum. 
Larger values of Bq tend to move this location downward; 
see Fig. HI A more detailed exploration of this is given by 
IKemeletal.l(l2012bl) . 

4.3.4. Comparison with the Parker instability 

NEMPI can be understood as a generalization of the Parker 
instability. This becomes evident when considering the stabil- 
ity criterion of NEMPI (RK07): 



(H p - H B ) 



dV, 



oil 



d/3 2 



>0, 



(30) 



/So 



where dV c s/df3 2 — |(1 — q p — dq p / 'din /3 2 ) and H B is the 
characteristic spatial scale of the mean magnetic field varia- 
tions. However, unlike the Parker instability, NEMPI can be 
excited even in a uniform mean magnetic field (Hb —t 00). 
The source of free energy of this instability is provided by the 
small-scale turbulent fluctuations. In contrast, the free energy 
in the Parker's magnetic buoyancy instability (Parker 1966) or 
in the interchange instability in plasma (Tserkovnikov 1960) 
is drawn from the gravitational field. In the absence of turbu- 
lence (q p = 0), condition d30b coincides with the criterion for 
the Parker's magnetic buoyancy instability (H p > Hb)- 

5. CONCLUSIONS 

Our DNS have shown that for an isothermal atmosphere 
with strong density stratification the total turbulent pressure 
is decreased due to the generation of magnetic fluctuations 
by the tangling of an imposed horizontal mean magnetic field 
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by the velocity fluctuations. This phenomenon strongly af- 
fects the mean Lorentz force so that the effective mean mag- 
netic pressure becomes negative. For our numerical model 
with approximately uniform turbulent rms velocity, the ratio 
of imposed to equipartition field strength changes with height, 
because the density decreases with height, while the imposed 
field is constant. This allows us to determine the full func- 
tional form of the effective mean magnetic pressure as a func- 
tion of normalized field strength for a single run. 

The form of the dependence of V c ff(B/B cq ) is similar to 
that found for simulations under rather different conditions 
(with or without stratification, with or without convection, 
etc). This dependence is found to be similar to that obtained 
earlier using both analytic theory (RK07) and direct numerical 
simulations (BKR), and the results are robust when changing 
the strength of the imposed field. 

In simulations where the turbulent velocity is nearly inde- 
pendent of height, the reduction of magnetic fluctuations oc- 
curs in the upper layers where the equipartition field strength 
decreases with height (models of type U). In models of type 
B, where the equipartition field strength is nearly constant 
in height, the magnetic fluctuations are found to be slightly 
stronger in the upper parts. 

In view of astrophysical applications, it is encouraging 
that <7 p o and /3 P seem to approach an asymptotic regime for 
R&M > 60. While it remains important to confirm this re- 
sult, a number of other aspects need to be clarified. Firstly, 
the issue of finite scale separation is important, i.e., the larger 
the wavenumber of structures in the mean field relative to kf, 
the less efficient the negative effective magnetic pressure will 
be. This needs to be quan tified. For example in the work of 
iBrandenburg et all d2011l) . where we had a scale separation 
ratio of 1:15, magnetic structures were best seen after averag- 
ing along the direction of the mean field. On the other hand, 
with a scale separation ratio of 1:30, st ructures where quit e 
pronounced already without averaging (Kernel e t aDl2012ah . 
In the Sun, the scale separation ratio is very large in the hori- 
zontal direction, but in the vertical direction the system is ex- 
tremely inhomogeneous and the vertical pressure scale height 
increases rapidly with depth. The significance of such effects 
on the negative effective magnetic pressure effect remains still 
quite unclear. 

Another aspect concerns the limitations imposed by the use 
of an isothermal equation of state. In the context of regular 
(non-turbulent) magnetic buoyancy, the system is known to 
be more unstable to the buoyancy instability when the fluctu- 
ations evolve isothermally (Acheson 1978; Hughes & Weiss 



11995b iMacGregor & Cassinellil 120031) . However, in the con- 
text of the negative effective magnetic pressure instability it is 
not clear in which direction this effect would work. The only 
simulations where the equation for specific entropy was taken 
into account are the simulations of Kapyla et al. (2012), who 
also considered an unstably stratified atmosphere. In their 
case, the negative effective magnetic pressure was found to 
be greatly enhanced (deeper minimum of V s s an d larger val- 
ues of B mm ). Again, this is a subject that deserves serious 
attention. 

The fact that the values of q p o and (3 p appear to be con- 
verged in the range 60 < Re^/ < 600 is significant, because 
this is also the regime in which small-scale dynamo action 
occurs. Small-scale dynamo action suppresses the magnetic 
pressure effect, which is the reason for the drop of q p o be- 
tween Re m of 40 and 60, but for larger ReM, the values of 
q p o seem roughly unchanged. 

In the present paper we have discussed applications mainly 
to the Sun. However, any hydromagnetic turbulence with 
strong stratification and large plasma beta should be sub- 
ject to the negative effective magnetic pressure phenomenon. 
Another relevant example might be accretion disks. Quasi- 
periodic oscillations and other light curve variations from ac- 
cretion disks have long been suspected to be caused by some 
kind of structures in thes e disks. Hydrodynamic vortices 
would be one possibility ( Abramowic z et alJ 119921). which 
could constitute long-lived structures (Barge & Sommeria 
11995b Hohansen et alJl2004t iLyra et alJl201 ll) . However, in 
view of the present results, structures caused by the nega- 
tive effective magnetic pressure instability might indeed be 
another candidate. 
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APPENDIX 
A. DERIVATION OF EQUATION (21) 



We use the quasi-linear appr oach or second order correlation approximation applied to a random flow with small magnetic and 
fluid Reynolds numbers (e.g. Moffatt 1978; Krause & Radler 1980; Rudigeret al. 2012). We eliminate the pressure term from 
the equation for the velocity fluctuations u by calculating V X (Vxit), rewrite the obtained equation a nd the induction equation 
for the magnetic fluctuations b in Fourier space, apply the two-scale approach (Roberts & Soward 1975), and neglect nonlinear 

terms, but retain molecular dissipative terms. This allows us to get the following equation for AIL from Eq. (0 in Fourier space: 



An£(fc,«) = - 



l\\ 



G* v {k,Lo) 
G v (k,u) 



nj (fe,w), 



(Al) 



where G u (k 7 uj) — {yk 2 + iu>) \ G ri (k, u) — (r]k 2 + iw)~ 
zero mean magnetic field) is given by 



L — Lr L ,Lr 7 jK Tn k n I3 r 



B n , and the background turbulence (with a 



n„' (fe,w) 



8irk 2 



(«. - *£) 






(A2) 
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Here E(k) is the energy spectrum, e.g., a power-law spectrum, E(k) = (q — l)(k/kf)~ q kj 1 with exponent 1 < q < 3 for 
the wavenumbers kf < k < kd, kf and kd are the forcing and dissipation wavenumbers, and we neglected for simplicity the 
anisotropy terms in Eq. (IA21 > which are proportional to A^ and XiXj, where A,; is a vector characterizing the anisotropy (see 
Appendix|B]l. We have taken into account that for small magnetic and hydrodynamic Reynolds numbers the small-scale dynamo 
is not excited, so that the background turbulence contains only the velocity fluctuations. We assume the frequency function $(w) 
to be a Lorentzian: <p(oj) = 1/[ttt c (u 2 + t~ 2 )]. This model for the frequency function corresponds to the correlation function 
(ui(t)uj(t + T)) oc exp(— r/r c ). After integration over u> and all angles in k space, and using Eq. (TT3l. we arrive at the following 
equations for q p0 and q s0 : 



Qpo 



QsO 



^2„,2 



15 



t;k 



1 + T r vk' 



k 2 E(k) 



15 



T c rjk 2 / (1 + t c vk 2 )(l + T c r\k 2 
k 2 E(k) 



■ dk, 



T c r)k 2 1 (1 + T c vk 2 )(l + T c r]k 2 ) 



dk, 



(A3) 



(A4) 



where we take into account that /3 <C 1. In the derivation of Eqs. ( IA3b and ( 1A41 i we used the following integrals for the integration 
in uj space: 



Crr? {-* v ^t ^r ^ 



(1 + T c vk 2 )(l + T c J]k 2 )' 



\Jff\ {-* ri {-* T {-* t QUI 



r]k 2 (l + r c r]k 2 ) ' 



(A5) 



We take into account that for small magnetic and fluid Reynolds numbers T c r)k 2 f 3> 1 and r c vk 2 - ^§> 1. 



where G T = (iu+r c 

In this limit the coefficients q p o and q s o, after integration over k, are given by: 



where 



where Prjv/ = Re^f /Re. 



Qpo 



C(q) 



15 



■Pr M )Rej\/Re, q s0 



dk 



15 



^ I T t 



(2 + Pr A . f )Re Af Re, 



1 - (k f /k d ) q " 



1 - (kf/kdY- 



(A6) 



(A7) 



B. THEORETICAL B DEPENDENCE OF q p , q s , AND q g 

In the following we summarize theoretical results for the B dependence of the coefficients q p , q s , and q g that enter in Eq. ( TT3b . 
These results were obtained for large magnetic and fluid Reynolds numbers using the r relaxation approach. We recall that the 
coefficient q p represents the isotropic turbulence contribution to the mean magnetic pressure, and q g is the anisotropic turbulence 
contribution to the mean magnetic pressure, while the coefficient q s is the turbulence contribution to the mean magnetic tension. 
We focus here on the case of anisotropic density-stratified background turbulence. Expressions for the isotropic case were given 
by RK07 and are summarized in BKR. Following BKR, we define (3 = B/B cq . We consider a plasma with a gas pressure that is 
much larger than the magnetic pressure, and the total pressure is always positive. 

We define the scale of the energy-carrying eddies as if w kj 1 . Due to density stratification, new terms emerge that are 

proportional to i 2 JH 2 These terms were absent in BKR, but otherwise the following formulae are identical. We also define the 

parameter e = (6q)/(mq), which takes into account the contributions caused by the small-scale dynamo (see RK07, where it was 
assumed for simplicity that the range of scales of magnetic fluctuations generated by the small-scale dynamo coincides with that 
of the velocity fluctuations). Table Q] suggests e = /3^ a ~ 0.3. 



For very weak mean magnetic fields, 4/3 <C Re A f , the values of 



<h(P) 



45 



(l + 91nRe M )(l-e) 



16 £ 2 
9flJ' 



8.03) 



2 P , q s , and q g are approximately constant and given by 

1 , , 8i 2 



15 



(l + 81nRe M )(l-e), q g (J3) 



3# P 2 ' 



-1/4 



for Re M ' « 4/3 C 1 we have 



ifi 16 £ 2 

q P W) = 25 [1 + 5| ln(4/3)| + 32 (3 2 \ (1 - e) + ^ 



1- 



16/3 2 



«.(/9) 



32 
15 



ln(4/3)| + - + 12/3 2 



(1-e), Qgtf) 



8i 2 

Jh 2 
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2, 



and for strong fields, 4/3 S> 1, we have 
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(B4) 
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Here we have taken into account that the anisotropic contributions to the nonlinear functions q p ((3) and q g (j3) for density-stratified 
background turbulence are given by 



Qg(P) 



■Qp(0) 



3 H I 



64/3 4 - A/3 2 + - + 



4/3 2 



fin 1 



arctan(Vo/3) 



8T2/3 3 



(B5) 



For the derivation of Eq. ( |B5t we used Eqs. (AlO)-(All) given by RK07 with the following model of the density-stratified 
background turbulence written in the Fourier space: 



(ui(k)uj(-k)) 



(u 2 )E(k) 



8nk 2 (k 2 + \ 2 ) 



Sij (k + A ) — ki kj — A, Xj + i (Aj kj ~ \j ki) 



where the velocity field satisfies the continuity equation in the anelastic approximation div u — m Aj, A, 
spectrum function is E(k) = (2/3) kj 1 (fc/fc/)" 5 / 3 for kf < k < kf Re 3/4 . 



(B6) 
Vip/p, the energy 



C. COMMENTS ON THE TEST-FIELD METHOD 



In the test- field method one uses a set of different test fields to determine all relevant components of the a and turbulent diffusiv- 
ity tensors. Furthermore, for finite scale separation ratios in space and time one also needs to represent all relevant wavenumbers 
and frequencies. The knowledge of all higher wavenumbers and frequenc ies allows one to compute th e integral kernels that de- 
scribe the nonlocality of turbulent transport; see Brandenburg et al. (2008c) for nonlocality in space and Hubbard & Brandenburg 
(2009) for nonlocality in time. The multitude of test fields does allow one to compute also those parts of the a and turbulent 
diffusivity tensors that do not enter in the particular problem at hand, but also those parts that enter under any other circumstances. 
An example is the evolution of a passive vector field where the same mean-field theory applies (Tilgner & Brandenburg 2008). 

Furthermore, given that we use the quasi-kinematic test-field method, we need to address the work of Courvoisier et al. (2010), 
who point out that this method fails if there is hydromagnetic background turbulence origin ating, for example, from small-s cale 
dynamo action. In such a case a fully nonlinear test-field method must be employed (see Rheinhardt & Brandenburg 2010, for 
details and implementation). However, it is worth noting that even in cases where small-scale dynamo action was expected, such 
as those of Brandenburg et al. (2008b) where values of Re a/ up to 600 were considered, the quasi-kinematic test- field method 
was still found to yield valid and self-consistent results, as was demonstrated by comparing the growth rate expected from the 
obtained coefficients of ag and rig. This gr owth rate was confirmed to be compatible with zero in the steady state. Finally, as 
shown in Rheinhardt & Brandenburg (2010), the quasi-kinematic method is valid if magnetic fluctuations result solely from an 
imposed field. In particular, the quasi-kinematic test-field method works even in cases in which magnetic fluctuations are caused 
by a magnetic buoyancy instability (Chatteriee et al. 2 01 lh . 

D. COMMENTS ON MEAN-FIELD BUOYANCY 



The work of Kitchatinov & Pipin ( 1993) is of interest in the pr esent context, becaus e it pre dicts the upward pumping of mean 
magnetic field. Here we discuss various aspects of this work. Kitchatinov & Pipin (1993) assumed that: (i) the gradient of 
the mean density is zero, (ii) the background turbulence is homogeneous, and (iii) the fluctuations of pressure, density and 
temperature are adiabatic. We also note that their analysis is restricted to low Mach number flows, although this is not critical for 
our present discussion. Since the gradient of the mean density is zero, the hydrostatic equilibrium, Vp = pg, exists only if the 
gradient of the mean temperature is not zero. This implies that the turbulent heat flux is not zero and temperature fluctuations are 

f;nerated by the tangling of this mean temperature gradient by the velocity fluctuations. Therefore, the key assumption made in 
itchatinov & Pipin ( 1993) that fluctuations of pressure, density and temperature are adiabatic, is problematic and the equation 
for the evolution of entropy fluctuations should be taken into account. This implies furthermore that the temperature fluctuations 
in Eq. (2.5) of their paper cannot be neglected. We avoid this here by considering flows with a non-zero mean density gradient 
and turbulence simulations that have strong density stratification. 
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